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Abstract 

We  consider  the  problem  of  tracking  multiple  maneuvering  targets  in  clutter  using  switching 
multiple  target  motion  models.  A  novel  suboptimal  filtering  algorithm  is  developed  by  applying  the 
basic  interacting  multiple  model  (IMM)  approach,  the  joint  probabilistic  data  association  (JPDA) 
technique  and  coupled  target  state  estimation  to  a  Markovian  switching  system.  The  algorithm  is 
illustrated  via  a  simulation  example  involving  tracking  of  two  highly  maneuvering,  at  times  closely 
spaced,  targets. 
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1  Introduction 


We  consider  the  problem  of  tracking  multiple  maneuvering  targets  in  clutter.  This  class  of  problem 
has  received  considerable  attention  in  the  literature  [2], [5], [6], [11].  The  switching  multiple  model 
approach  has  been  found  to  be  quite  effective  in  modeling  highly  maneuvering  targets  [2],  [4]- [7]. 
In  this  approach  various  “modes”  of  target  motion  are  represented  by  distinct  kinematic  models, 
and  in  a  Bayesian  framework,  the  target  maneuvers  are  modeled  by  switchings  among  these  models 
controlled  by  a  Markov  chain.  In  the  presence  of  clutter,  the  measurements  at  the  sensors  may  not 
all  have  originated  from  the  target-of-interest.  In  this  case  one  has  to  solve  the  problem  of  data 
association.  An  effective  approach  in  a  Bayesian  framework  is  that  of  probabilistic  data  association 
(PDA)  [2],  [4]  for  a  single  target  in  clutter  and  that  of  joint  probabilistic  data  association  (JPDA) 
[2],  [5],  [11]  for  multiple  targets  in  clutter. 

It  is  assumed  that  the  number  of  targets  is  known  (say  N)  and  that  for  each  target,  a  track  has 
been  formed  (initiated)  and  our  objective  is  that  of  track  maintenance.  In  [12]  such  a  problem  has 
been  considered  for  a  single  target  using  multiple  sensors,  PDA  and  switching  multiple  models.  The 
optimal  solution  (in  the  minimum  mean-square  error  sense)  to  target  state  estimation  given  sensor 
measurements  and  absence  of  clutter,  requires  exponentially  increasing  (with  time)  computational 
complexity;  therefore,  one  has  to  resort  to  suboptimal  approximations.  For  the  switching  multiple 
model  approach,  the  interacting  multiple  model  (IMM)  algorithm  of  [7]  has  been  found  to  offer  a 
good  compromise  between  the  computational  and  storage  requirements  and  estimation  accuracy 
[13].  In  the  presence  of  clutter,  one  has  to  account  for  measurements  of  uncertain  origin  (target 
or  clutter?).  Here  too,  in  a  Bayesian  framework,  one  has  to  resort  to  approximations  to  reduce 
the  computational  complexity,  resulting  in  the  PDA  filter  [2],  [5],  [1],  [12].  In  [12]  the  IMM  algo¬ 
rithm  has  been  combined  with  the  PDA  filter  in  a  multiple  sensor  scenario  to  propose  a  combined 
IMM/MSPDAF  (interacting  multiple  model/  multisensor  probabilistic  data  association  filter)  algo¬ 
rithm.  In  [2]  and  [11]  multiple  targets  in  clutter  (but  without  using  switching  multiple  models)  have 
been  considered  using  JPDA  filter  which,  unlike  the  PDA  filter,  accounts  for  the  interference  from 
other  targets.  Various  versions  of  IMMJPDA  filters  for  multiple  target  tracking  using  switching 
multiple  models  may  be  found  in  [3],  Sec.  6  of  [5],  [9]  and  [10].  While  [9]  and  [10]  present  uncoupled 
filters  (i.e.  assume  that  different  target  states  are  mutually  independent  conditioned  on  the  past 
measurements),  [3]  and  [5]  present  IMMJPDA  coupled  filters  where  the  conditional  target  state 
independence  assumption  is  not  made.  It  has  been  noted  in  {8]  that  the  IMMJPDA  coupled  filter 
equations  of  [3]  and  [5]  are  heuristic.  [8]  presents  an  “exact”  JPDA  coupled  filter  for  non-switching 
models  using  the  framework  of  a  linear  descriptor  system. 

In  this  correspondence  we  extend  the  approach  of  [9]  which  pertains  to  uncoupled  filtering,  to 
IMMJPDA  coupled  filtering.  The  only  approximations  made  are  those  typical  for  IMM  approaches; 
there  are  no  other  heuristics  as  in  [3]  and  [5] .  We  use  the  standard  Markovian  switching  state-space 
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systems  which,  as  discussed  in  Remark  4  in  Sec.  3,  is  equivalent  to  the  linear  descriptor  system 
framework  of  [8],  Track  initialization  (formation)  is  assumed  to  have  been  made  for  each  target. 
“Standard”  assumptions  are  used  for  JPDAF  ([11],  p.  310  of  [5]):  a  measurement  can  have  only 
one  source;  among  the  possibly  several  validated  measurements,  at  most  one  of  them  can  be  target- 
originated  and  the  remaining  validated  measurements  are  assumed  to  be  due  to  false  alarms  or 
clutter,  and  are  modeled  as  independently  and  identically  distributed  (i.i.d.)  with  uniform  spatial 
distribution  over  the  entire  validation  region  (“across  all  targets”).  As  in  [12]  and  [9]  we  use 
sequential  updating  of  the  state  estimates  with  sensor  measurements.  As  noted  in  Remark  1  in 
Sec.  3,  this  is  a  suboptimaJ  approach  since  joint  measurement  association  events  across  sensors  are 
not  taken  into  account;  for  an  optimal  approach,  one  should  follow  [14]-[16]  -  see  Remark  1  for 
further  details. 

2  Problem  Formulation 

Assume  that  there  are  total  N  targets  with  the  target  set  denoted  as  Tjv  :=  {1, 2,  •  •  • ,  N}.  Assume 
that  the  dynamics  of  each  target  can  be  modeled  as  one  of  the  n  hypothesized  models.  The  model 
set  is  denoted  as  Mn  :=  {1, 2,  •  ■  •  ,n}  and  there  are  total  q  sensors.  For  target  r  (r  G  TJv),  the  event 
that  model  i  is  in  effect  during  the  sampling  period  (tfc_i,tjfc]  will  be  denoted  by  M^(r),  Although 
all  the  targets  share  a  common  model  set,  any  two  targets  may  be  in  different  motion  status  from 
time  to  time.  For  the  j-th  hypothesized  model  (mode),  the  state  dynamics  and  measurements  of 
target  r  (r  G  TJv)  are  modeled  as 

^kir)  =  Fl_i{r)xk-iir) -h  (1) 

zi(r)  =  (xk  (r))  +  (r)  for  /  =  1 ,  •  •  • ,  q  (2) 

where  Xk{r)  is  the  system  state  of  target  r  at  and  of  dimension  rix  (assuming  all  targets  share 
a  common  state  space),  z\.{r)  is  the  (true)  measurement  vector  (i.e.  due  to  target  r)  from  sensor 
I  at  tk  and  of  dimension  n^/,  and  Gi_j(r)  are  the  system  matrices  when  model  j  is  in 

effect  over  the  sampling  period  (tfe-i,tfc]  for  target  r  and  is  the  nonlinear  transformation  of 
Xk(r)  to  zl(r)  (I  =  1, 2,  •  •  • ,  g)  for  model  j.  The  process  noise  (r)  and  the  measurement  noise 
(r)  are  mutually  uncorrelated  zero-mean  white  Gaussian  processes  with  covariance  matrices 
QI_i  (same  for  all  targets)  and  (same  for  all  targets),  respectively.  At  the  initial  time  to,  the 
initial  conditions  for  the  system  state  of  target  r  under  each  model  j  are  assumed  to  be  Gaussian 
random  variables  with  the  known  mean  XQ(r)  and  the  known  covariance  Po(r).  The  probability  of 
target  r  in  model  j  at  to,  fJ-oir)  =  P{M^{r)},  is  also  assumed  to  be  known.  The  switching  from 
model  A/fe_i(r)  to  model  M^(r)  is  governed  by  a  finite-state  stationary  Markov  chain  (same  for 
all  targets)  with  known  transition  probabilities  pij  =  P{Mi{r)\Ml_■^^{r)}.  Henceforth,  tk  will  be 
denoted  by  A:. 
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In  coupled  state  estimation  the  states  of  all  targets  are  estimated  jointly  [5].  To  this  end  define 
the  “global”  state 

Xfc  :=  col{xjfe(l),  Xfc(2),  ••  •  ,a;fc(Ar)}  (3) 

and  the  corresponding  matrices/vectors  J  :=  col{ji,  ja,  •  •  •  ,jN}  where  €  Mn  is  model  j  for  tar¬ 
get  m,F/  :=  block -diag{i^>(l),  ;=  block  -  diag  {cini),  G^(2),  •••, 

v(  :=  col  Vfc^(2),  •  •  •  Then  we  have  the  state  equation  for  the  N  tar¬ 

gets  as 

=  ■f’fc-i  Xk-i  +  G(_  I  vl_  1  (4) 

where  E{vlvl'}  =  Ql  .=  block  -  diag  •  ••  Similarly  define  the  global  measurement 

vector  at  sensor  1  as 


zi:=col{4(l),4(2),...,4(Ar)}  (5) 

and  the  corresponding  vectors  h^'^{xk)  '.=  col 

:=  col  ^k(2),  ■  •  •  where  E{w^^w^^'}  =  :=  block  -  diag 

Then  the  measurement  equation  for  N  targets  at  sensor  I  (assuming  no  clutter  and  perfect  detec¬ 
tions)  is  given  by 

zl.  =  h‘^’^{xk)+wi'^  for  f  =  ,  (6) 


Define  the  global  mode  :=  •  •  • ,  mI^{N)}.  The  various  targets  are  assumed  to  evolve 

independently  of  each  other.  Therefore,  the  transition  probability  for  the  global  modes  are  given 
by 

J=i 

Similarly  we  have 

:=  P{mA(1),  Af«(iV))  =n,4'(i).  (8) 

/=1 


Regarding  the  measurements  at  sensor  I,  we  follow  the  notations  and  definitions  used  used 
in  [9].  At  any  time  k,  some  measurements  may  be  due  to  clutter  and  some  due  to  the  target. 
The  measurement  set  (not  yet  validated)  generated  by  sensor  /  at  time  k  is  denoted  as  Z^.  := 

where  mi  is  the  number  of  measurements  generated  by  sensor  I  at  time 
k.  Variable  (i  =  1,  is  the  ith  measurement  within  the  set.  The  validated  set  of 

measurements  of  sensor  I  at  time  k  will  be  denoted  by  Y^,  containing  mi  (<  mi)  measurement 
vectors.  The  cumulative  set  of  validated  measurements  from  sensor  I  up  to  time  k  is  denoted  as 
‘  The  cumulative  set  of  validated  measurements  from  all  sensors  up  to  time 
k  is  denoted  as  2*  =  •  •  • , where  q  is  the  number  of  sensors. 
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Assuming  there  are  no  unresolved  measurements  (i.e.  measurement  associated  with  two  or  more 
targets  simultaneously),  any  measurement  therefore  is  either  associated  with  a  target  or  caused  by 
clutter.  Our  goal  is  to  find  the  global  state  estimate  '•=  E{xk\2i}  and  the  associated  error 
covariance  matrix  Pk\k  =  E{[xk  —  where  x'k  denotes  the  transpose  of  Xk- 

Included  in  the  above  formulation  is  state  estimates  of  individual  targets. 


3  IMM/JPDA  Coupled  Filtering  Algorithm 


We  now  modify  the  IMM/JPDA  algorithm  of  [9]  to  apply  to  the  coupled  system  (3)-(8);  it  will  be 
called  IMM/JPDACF  (CF  stands  for  coupled  filter).  The  approach  of  [9],  in  turn,  is  based  on  the 
approaches  of  [12],  [11],  [5]  and  [2].  As  in  [12]  and  [9],  for  convenience,  we  confine  our  attention 
to  the  case  of  2  sensors;  however,  the  algorithm  can  be  easily  adapted  to  the  case  of  arbitrary  q 
sensors.  As  the  IMM/MSPDAF  algorithm  is  well-explained  in  [12]  and  Sec.  4.5  of  [5],  the  JPDAF 
algorithm  is  well-explained  in  Sec.  6.2  of  [5]  and  Sec.  9.3  of  [2]  and  the  IMM/JPDA  filter  is  given  in 
detail  in  [9],  we  will  only  briefly  outline  the  basic  steps  in  “one  cycle”  of  the  IMM/JPDA  coupled 
filter. 

Assumed  available:  Given  the  state  estimate  :=  E{xk-i\M^_^,Z\~^),  the  associated 

covariance  Pk-\\k-\  conditional  mode  probability  =  P[M/_j  [2*“^]  at  time  fc  —  1  for 

each  global  mode  J  £  .Ad„  :=  Mn  x  •  •  •  x  A4„. 

Step  3.1.  Interaction  —  mixing  of  the  estimate  from  the  previous  time  (VJ  e  Mn)’- 


predicted  mode  probability:  P{Mk\Z^  (9) 

/ 

mixing  probability:  (10) 

mixed  estimate:  =  (H) 

I 

covariance  E{[xk-\  -  2fc£i|fe_i][2:fe-i  -  ^-i\k~i]'Wk  ^ ^i~^}  of  fh®  mixed  estimate: 

/ 

Step  3.2.  Predicted  state  and  measurements  for  Sensor  1  (VJ  €  Mn)- 

State  prediction:  :=  E{xk\M^ ,Z\-^}  =  Fi_-^^k-^\k-v  (13) 


State  prediction  error  covariance  E{[xk  —  xLk--^[^k  — 


pJ 

^k\k-l 


nJ  pOJ  ipJ'  I  /nf J  /nfj' 
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Using  (2)  and  (13),  the  global  mode-conditioned  predicted  measurement  for  sensor  1  is  := 
Using  linearized  (6),  the  covariance  of  the  mode-conditioned  residual 
:=  col{2j^*‘\  •  •  • ,  is  given  by 


:=  Zu'  - 


where  is  the  first  order  derivative  (Jacobian  matrix)  of  at  Note  that  (15) 

assumes  that  originates  from  the  target  r. 

Step  3.3.  Measurement  validation  for  sensor  1  (Vj  €  Mn)-  First  perform  measurement 
validation  for  each  target  r  (r  €  Tn)  separately.  For  target  r,  the  validation  region  is  taken  to  be 
the  same  for  all  models,  i.e.,  as  the  largest  of  them.  Let  Sj!^{r)  denote  the  Uzi  x  rizi  submatrix  of 
including  the  rows  and  columns  of  the  latter  numbered  as  (r  —  l)nzi  -l-m,  m  =  1, 2,  •  •  • ,  r.  That 
is,  iS'^’^(r)  is  based  on  the  information  relevant  to  target  r  only.  Let  ^’’’^(r)  denote  the  Uzi  x  1  sub¬ 
column  of  z^’^  including  the  rows  of  the  latter  numbered  as  (r  —  l)nzi  +m,  m  =  1,2,  -  •  •  ,r;  that  is, 
^"■’^(r)  is  the  mode-conditioned  predicted  measmrement  of  target  r  for  sensor  1.  Let  (|A|  =  det(.(4)) 

3r  :=  arg  |  m^  [5^’^  (r)  |  J  .  (16) 


^  (  i  =  1, 2,  •  •  • ,  mi)  is  validated  if  and  only  if 


Then  measurement  ^  — 


(17) 


where  7  is  an  appropriate  threshold.  The  volume  of  the  validation  region  with  the  threshold 
7  is  V^{r)  :=  where  Uzi  is  the  dimension  of  the  measurement  and  is 

the  volume  of  the  unit  hypersphere  of  this  dimension.  After  performing  the  validation  for  each 
target  separately,  the  volume  of  validation  region  for  the  whole  target  set  is  approximated  by 
=  E^i 

Step  3.4.  State  estimation  with  validated  measurements  from  sensor  1  (VJ  6  Mn)' 
Given  Zl  :=  ‘  define  the  set  of  validated  measurement  for  sensor  1  at  time 

k  as  :=  Vk^'^^y ' '  ‘  where  mj  is  total  number  of  validated  measurement  for  sensor 

1  at  time  k.  and  :=  with  1  <  /i  <  Z2  <  •  •  •  <  Imi  <  when  mi  ^  0.  We  now  consider 
joint  probabilistic  data  association  across  targets  following  [5]  and  [2],  but  for  global  target  state; 
note  that  here  we  consider  only  sensor  1  (see  also  Remark  1  later).  A  marginal  association  event 
Bir  is  said  to  be  effective  at  time  k  when  the  validated  measurement  is  associated  with  (i.e. 
originates  from)  target  r  (r  =  0,  -  ■  •  ,N  where  r  —  0  means  that  the  measurement  is  caused  by 
clutter).  Assuming  that  there  are  no  unresolved  measurements,  a  joint  association  event  ©  is 
effective  when  a  set  of  marginal  events  {Bir}  holds  true  simultaneously.  That  is,  ©  =  ^^lBir^ 
where  rj  is  the  index  of  the  target  to  which  measurement  is  associated  in  the  event  under 
consideration.  Define  the  validation  matrix  (as  in  [5]  and  [2]) 

D=[wir]  i  =  !,••■, mi,  r  =  0,’--,N  (18) 
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where  Wjr  =1  if  the  measurement  i  lies  in  the  validation  gate  of  target  r,  else  it  is  zero.  A  joint 
association  event  ©  is  represented  by  the  event  matrix 

f2(©)  =  [uiir(©)]  r  =  0,---,iV  (19) 

where  a)ir(©)  =  1  if  9ir  C  ©,  =0  otherwise.  A  feasible  association  event  is  one  where  a  measurement 
can  have  only  one  source  ^ir(©)  =  1  Vi)  and  where  at  most  one  measurement  can  originate 
from  a  target.  The  feasible  association  joint  events  ©  are  mutually  exclusive  and  exhaustive. 

Following  the  definitions  in  [5]  and  [2],  define  the  binary  measurement  association  indicator 
ri(©)  :=  S^jiDirC©),  (i  =  1,  •••,7ni),  to  indicate  whether  the  validated  measurement  is 
associated  with  a  target  in  event  ©.  Further,  the  number  of  false  (unassociated)  measurements  in 
event  ©  is  0(©)  =  We  will  limit  our  discussion  to  nonparametric  JPDA  [5]  and 

[2].  One  can  evaluate  the  hkelihood  that  the  global  mode  is  J  as 

e 

(20) 

e 

where  (Sec.  6.2  of  [5],  Sec.  9.3  of  [2]) 

=  <5.(©):=f]a)i,(©),  (21) 

*=1  i=0 

Pp  is  the  detection  probability  at  sensor  1  (assumed  to  be  the  same  for  all  targets)  and  e  >  0  is 
a  “diflFuse”  prior  (for  nonparametric  modeling  of  clutter)  whose  exact  value  is  irrelevant.  Unlike 
[9] ,  we  do  not  assume  that  the  states  of  the  targets  (including  the  modes)  conditioned  on  the  past 
observations  are  mutually  independent.  Then  we  have 

where  i^^(©)  C  subset  of  the  validated  measurements  yj|.\  consisting  of  the  measurements 

associated  with  the  targets  as  specified  by  ©.  The  number  of  measurements  in  Y^  (©)  equal  mi  — 
<f>{Q)  where  <f>{Q)  is  the  number  of  false  alarms. 

Define  a  mi  x  [mi  —  <^(©)]  matrix  £!(©)  as  a  submatrix  of  fl(©)  obtained  by  deleting  the  first 
column  and  all  null  columns  of  Q{@).  Then  for  a  given  ©,  we  have  a  measurement  vector  1)^  (©) 
of  dimension  (X)S:i  'n(©))^2i  given  by 

n'(©)  =  ® 0'(e))col{t/f  \  r  =  1, 2, •  •  •  ,mi}  (23) 

where  we  stack  up  all  target-associated  validated  measurements  in  ©  in  ascending  order  of  targets, 
In  is  the  n  X  n  identity  matrix,  and  the  symbol  0  denotes  the  Kronecker  product.  Define  a 
[(mi  —  0(©))7i2i]  X  [N^x]  matrix  as  a  submatrix  of  obtained  by  deleting  all  i-th  block 
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rows  {rizi  X  .)  of  for  which  5i(©)  =  0.  That  is,  we  have  modified  to  keep  only  the  block 
elements  associated  with  target-associated  measurements  in  0.  It  then  follows  that  the  linearized 
measurement  equation  for  (0)  is  given  by 

=  (24) 


Conditioned  on  the  joint  association  event  0  and  mode  J,  the  “coupled”  innovations  is  given 


by 


FfcH©)  -  (©)  if  =  1  for  some  r  e  {1, 2, •  •  • , N), 

0  otherwise, 


(25) 


where  z^'^(G)  is  a  subvector  of  obtained  by  deleting  all  i-th  block  rows  {uzi  x  1)  of  for 
which  5t(0)  =  0.  The  conditional  pdf  (probability  density  function)  of  the  validated  measurements 
(0)  given  their  origins  (specified  by  0)  and  the  global  target  mode  J,  is  given  by 


p[yfci(0)|M/,2i"-i]  =  Ar(ni(©);2^’i(©),5f  (0))  (26) 

where 

■V(i;!/,P)  ;=  |2irPr‘''^exp|-i(x-  i/)'P"‘(x  -y)]  (27) 

S"(e)  ;=  (0)P4-,if^'{e)'  +  (28) 


The  probability  of  the  joint  association  event  0  given  that  global  mode  J  is  effective  from  time 
k  —  1  through  fc  is 

Pi\Q)-=p{Q\Mk^zt\yk}  =  \p\yk\^M.zt-^]pm  (29) 

where  the  first  term  can  be  calculated  from  (22)-(28),  the  second  term  from  (21),  and  c  is  a 
normalization  constant  such  that  P{^\^k  •> 

Using  (from  (13))  and  its  covariance  (from  (14)),  one  computes  the  partial  update 

and  its  covariance  P^l  following  the  standard  PDAF  [5],  [2],  except  that  the  global  state  is 
conditioned  on  0,  not  the  marginal  events  Oir',  details  follow. 

Kalman  gain:  H'/(e)  = (e)'|S«(e)]->. 

Partial  update  of  the  state  estimate:  ®^]fc(©)  :=  E{xk\Q, ,Y^} 

^k\k-i  "b  (®)^fc’^(®)  ^f  <^r(®)  7^  9  for  some  1  <  r  <  iV 
H\k-i  if  <!.(©)  =  0Vr€  {1,2,...,  AT}. 

:=  (e)x"(e) 

e 
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(30) 

(31) 

(32) 


Covariance  of  Sfii  :  =  JiV,  -  Y,  ft"'‘(e)W'/(e)sf 

©:©/©o 

+  '£Pi\e)Wi^{e)ui’He)ui'\ey  Wilier 

© 


j;/3f(e)»'j'(e)of(e) 
.  © 


.  © 


(33) 


where  ©o  denotes  ©  for  which  (Jr(©)  =  0  Vr  e  (1, 2,  •  •  • ,  N}.  Eqn.  (33)  follows  in  a  manner  similar 
to  eqn.  (3.4.2-10)  in  [5]. 


Step  3.5.  The  mode-conditioned  predicted  measurements  for  sensor  2  (VJ  G  Mn)'-  The 
“predicted”  measurement  for  sensor  2  is  given  by 


(34) 

The  covariance  of  the  global  mode-conditioned  residual  ~  ^  given  by 

where  is  the  first  order  derivative  (Jacobian  matrix)  of  h‘^’^(.)  at  . 

Step  3.6.  Measurement  validation  for  sensor  2;  This  is  similar  to  Step  3.3  where  we  replace 
with  with  mi  with  m2,  Vj^ir)  with  V^{r),  and  1^^  with  V^.  Details  are  similar 

to  that  in  Step  3.3,  hence  are  omitted. 


Step  3.7.  Update  with  validated  measurements  for  sensor  2  (VJ  6  Mn)'-  This  step  is 
similar  to  Step  3.4.  Using  the  validated  measurements  obtained  from  Step  3.6  and  starting  from 
and  one  computes  the  final  updates  and  and  the  likelihood 

(se) 


The  details  are  similar  to  that  in  Step  3.4,  hence  are  omitted. 

Step  3.8.  Update  of  mode  probabilities  (Vj  G  Mn,  Vr  G  Tn): 

4:=PlMl\Zi]=^4-Ai‘Af  (37) 

where  c  is  a  normalization  constant  such  that  Y,j  t^k  =  For  individual  targets  we  have 

Mi»:=rK'(r)|Zfl=  f;...  ^  ^  (3g) 

jl=l  Jr-l  =  ljr+l=l  iiV  =  l 

with  J  =  (ji,  •  •  • ,  jjv)  in  (37). 

Step  3.9.  Combination  of  the  mode-conditioned  estimates  (Vr  G  T^):  The  final  global 
state  estimate  update  at  time  k  is  given  by  =  Yi,j^k\kl^k  its  covariance  is  given  by 
Pk\k  =  T,j{Pk\k  +  [H\k~^k\k][H\k~^k\k]']  ^^k■  The  state  estimate  £fc|fc(r)  for  target  r  is  the 
nx-subvector  of  %|i.  consisting  of  elements  (r  —  l)nx  +  m,  m  =  1, 2,  •  •  • ,  Ux. 
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Remark  1.  In  the  above  algorithm  we  used  sequential  updating  of  the  state  estimates  with 
measurements  (one  sensor  at  a  time  -  see  Steps  3.4  and  3.7)  as  in  [9]  and  [12].  This  approach  is 
suboptimal  but  leads  to  computational  savings  as  one  does  not  have  to  simultaneously  associate 
measurements  across  sensors  (as  in  [14]-[16]).  In  Step  3.4  we  are  interested  in  (an  approximation 
to)  E{xk\M^,  which  is  decomposed  as  in  (32)  conditioned  on  ©’s;  measurements  are 

not  considered  in  this  step.  If  one  were  to  seek  E{xk\M^ , ,Y^ ,Y^},  then  we  would  have  to 
follow  the  approach  of  [14]- [16]  by  picking  all  possible  association  events  across  sensors  also. 

Remark  2.  Compared  to  the  uncoupled  filtering  of  [9]  where  the  equations  are  formulated 
conditioned  on  marginal  association  events  6ir,  here  we  have  conditioning  on  joint  association  events 
©  for  coupled  filtering.  Eqn.  (26)  does  not  decompose  into  the  product  of  marginal  probabilities 
as  in  [9]. 

Remark  3.  Partition  the  set  of  all  ©s  into  disjoint  sets  ©jS  such  that 
:=  {0 1  J^(©)  =  5^(©)  Vr,  ©  e  ©i} 

where  i  =  -  ,K.  For  instance,  for  =  2,  we  have  K  =  4  with  ©i  =  all  ©s  in  which  there 

are  two  validated  measurements  associated  with  two  targets,  ©2  =  all  ©s  in  which  one  validated 
measurement  is  associated  with  target  1  and  none  with  target  2,  ©3  =  all  ©s  in  which  one  validated 
measurement  is  associated  with  target  2  and  none  with  target  3,  and  ©4  =  all  ©s  in  which  none 
of  the  validated  measurements  are  associated  with  any  target.  It  is  then  easily  seen  that 

S^’^{Q)  and  /3^’^(©)  in  Step  3.4,  all  are  invariant  for  ©  €  ©i.  This  fact  can  be  used  to 
simplify  computations  in  (31)-(33).  Similar  comments  apply  to  Step  3.7. 

Remark  4.  If  one  substitutes  (23)  into  (24),  then  one  obtains  a  linear  descriptor  system  type 
of  equation  such  as  (12)  in  [8].  Therefore,  the  standard  state-space  system  framework  used  in  this 
paper  and  the  linear  descriptor  system  framework  used  in  [8]  are  equivalent  (except  that  [8]  uses 
non-switching  models  whereas  we  use  Markovian  switching  models).  We  have  retained  the  notation 
and  formulation  of  the  earlier  papers  in  the  field  (e.g.  [5],  [7],  [9]  and  [12])  whereas  [8]  prefers  to 
follow  a  linear  descriptor  system  formulation. 

4  Simulation  Examples 

We  now  consider  tracking  two  highly  maneuvering  targets  in  clutter. 

4.1  Example  1 

The  True  Trajectory:  Target  1  starts  at  location  [21689  IO84O  40]  in  Cartesian  coordinates  in 
meters.  The  initial  velocity  (in  m/s)  is  [—8.3  —  399.9  0]  and  the  target  stays  at  constant  altitude 
with  a  constant  speed  of  400m/s.  Its  trajectory  is:  a  straight  line  with  constant  velocity  between 
0  and  17s,  a  coordinated  turn  (0.15  rad/s)  with  constant  acceleration  of  60  m/s^  between  17  aiid 
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30s,  a  straight  line  with  constant  velocity  between  30  and  55s,  a  coordinated  turn  (0.1  rad/s)  with 
constant  acceleration  of  40  m/s^  between  55  and  70s,  and  a  straight  line  with  constant  velocity 
between  70  and  87s.  Target  2  starts  at  location  [30000  -  3040  40]  in  Cartesian  coordinates  in 
meters.  The  initial  velocity  (in  m/s)  is  [—382  157  0]  and  the  target  stays  at  constant  altitude  with 
a  constant  speed  of  413m/s.  Its  trajectory  is:  a  straight  line  with  constant  velocity  between  0  and 
44s,  a  coordinated  turn  (0.075  rad/s)  with  constant  acceleration  of  30  m/s^  between  44  and  59s, 
and  a  straight  line  constant  velocity  between  59  and  87s. 

The  Target  Motion  Models;  These  are  exactly  as  in  [9].  The  motion  models  for  the  two 
targets  are  identical.  In  each  mode  the  target  dynamics  axe  modeled  in  Cartesian  coordinates 
as  Xk{r)  =  F{r)xic-i{r)  +  G{r)vk-i{r)  where  the  state  of  the  target  is  position,  velocity  and 
acceleration  in  each  of  the  3  Cartesian  coordinates  (r,  y  and  z).  Model  1  for  nearly  constant 
velocity  model  with  zero  mean  perturbation  in  acceleration;  Model  2  for  Wiener  process  acceleration 
(nearly  constant  acceleration  motion);  Model  3  for  Wiener  process  acceleration  (model  with  large 
acceleration  increments,  for  the  onset  and  termination  of  maneuvers).  The  details  regarding  these 
models  may  be  found  in  [9].  The  initial  model  probabilities  for  two  targets  axe  identical:  =  0.8, 

fXQ  =  0.1  and  ^0  =  0.1.  The  mode  switching  probability  matrix  for  two  targets  is  also  identical  and 
is  as  in  [9]. 

The  Sensors:  Two  sensors  (we  assume  collocation,  and  time  synchronization  of  observations, 
etc.)  are  used  to  obtain  the  measurements.  The  measurements  from  sensor  I  for  model  j  are 
z]^  =  h^’'’{xk)+'u^k  ,  Z  =  1, 2,  reflecting  range  and  azimuth  angle  for  sensor  1  (radar),  and  azimuth  and 
elevation  angles  for  sensor  2  (infrared).  The  range,  azimuth  and  elevation  transformations,  respec¬ 
tively,  are  given  by  r  =  (r^  +y^  +  a  =  tan~^(y/x),  e  =  lan~^[z/{x'^  -f  The  mea^ 

surement  noise  for  sensor  I  is  assumed  to  be  zero-mean  white  Gaussian  with  known  covariances 
E}  =  diag[g'r,9ai]  =  diag[400m^,  49mrad^]  with  qai  and  qr  denoting  the  variances  for  the  radar  az¬ 
imuth  and  range  measurement  noises,  respectively,  and  S?  =  diag[go2)9e]  =  diag[4mrad^,  4mrad^] 
with  qa2  and  qe  denoting  the  variances  for  the  infrared  sensor  azimuth  and  elevation  measurement 
noises,  respectively.  Both  sensors  are  assumed  to  be  located  at  the  coordinate  system  origin.  The 
sampling  interval  was  T  =  Is  and  it  was  assumed  that  the  probability  of  detection  Pd  =  0.997  for 
both  sensors. 

The  Clutter;  For  generating  false  measurements  in  simulations,  the  clutter  was  assumed  to 
be  Poisson  distributed  with  expected  number  of  Ai  =  20  x  10“®/m  mrad  for  sensor  1  and  A2  = 
2  X  10“'^/mrad^  for  sensor  2. 

Other  Parameters:  The  gates  for  setting  up  the  validation  regions  for  both  the  sensors  were 
based  on  the  threshold  7  =  16  corresponding  to  a  gate  probability  Pq  =  0.9997. 

Simulation  Results:  The  results  were  obtained  from  1000  Monte  Carlo  runs.  Fig.  1  shows  the 
true  trajectory  of  the  two  targets  and  the  distance  between  the  two  targets  as  a  function  of  time. 
The  two  targets  start  out  far  apart,  move  close  to  each  other  from  30  to  45  seconds,  and  then 
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move  apart  again.  Fig.  2  shows  the  results  of  the  proposed  IMM/ JPDACF  based  on  993  successful 
runs  (target  swapping  occurred  in  7  out  of  1000  runs).  Fig.  3  shows  the  results  of  the  uncoupled 
IMM/JPDAF  of  [9]  based  on  982  successful  runs  (target  swapping  occurred  in  18  out  of  1000 
runs).  Comparing  Figs.  2  and  3  we  see  that  when  there  is  no  target  swapping,  differences  between 
the  coupled  and  uncoupled  filters  are  small.  To  assess  the  computational  requirements  of  the  two 
approaches,  we  computed  the  CPU  time  needed  to  execute  one  time  step  (averaged  over  10  runs 
and  87  time  steps  in  each  run)  in  MATLAB  6.5  on  a  1  GHz  (Mobile)  Pentium  III  operating  under 
Windows  XP  (professional).  The  uncoupled  IMM/JPDAF  needs  0.1063  secs  compared  to  0.1369 
secs  required  by  IMM/JPDACF.  Thus  with  a  29%  increase  in  computational  cost,  IMM/JPDACF 
results  in  a  61%  fewer  target  swappings  compared  with  uncoupled  IMM/JPDAF.  Finally,  neither 
approach  experienced  any  loss  of  tracks,  only  target  swapping. 

4.2  Example  2 

We  now  consider  the  same  secnario  as  that  in  Example  1  except  for  a  linear  shift  in  the  y-direction 
in  the  trajectory  of  target  2.  Target  2  now  starts  at  [30000  —  3040  +  d  40]m  with  d=  -500,  -250, 
250  or  500m.  When  d  =  0,  we  get  Example  1.  Different  values  of  d  lead  to  different  separations 
and  encounters  between  the  trajectories  of  the  two  targets:  Fig.  5(a), (b)  shows  the  true  trajectory 
of  the  two  targets  for  two  different  values  of  d  (-500m  and  500m).  For  each  value  of  d,  results 
were  obtained  from  100  Monte  Carlo  runs.  Table  1  shows  the  number  of  successful  runs  (no  target 
swapping)  versus  d  for  the  two  approaches  IMM/JPDACF  and  IMM/JPDAF.  Fig.  4(c)-(f)  shows 
the  position  error  versus  time  for  different  values  of  d  averaged  over  only  the  successful  runs.  It  is 
seen  from  Table  1  and  Fig.  4  that  IMM/JPDACF  is  either  better  than  IMM/JPDAF  (e.g.  d=-250m 
or  Om)  or  similar  to  IMM/JPDAF  (other  values  of  d:  larger  separation)  in  performance. 

5  Conclusions 

We  proposed  a  novel  IMM/JPDA  coupled  filtering  algorithm  for  state  (position,  velocity  and  ac¬ 
celeration)  estimation  for  multiple  highly  maneuvering  targets  in  clutter.  The  algorithm  was  illus¬ 
trated  via  a  simulation  example  where  with  a  29%  increase  in  computational  cost,  IMM/JPDACF 
resulted  in  a  61%  fewer  target  swappings  compared  with  uncoupled  IMM/JPDAF;  neither  approach 
experienced  any  loss  of  tracks. 
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The  trajectory  of  the  two  teraets 


Figure  1:  Example  1.  The  true  trajectories  of  the  maneuvering  targets  (read  left-to-right,  top- 
to-bottom):  (a)  Position  in  the  xy— plane,  (b)  x  and  y  velocities,  (c)  x  and  y  accelerations,  (d) 
distance  between  the  targets. 


d 

IMM/JPDACF 

IMM/JPDAF 

-500m 

100/100 

100/100 

-250m 

99/100 

96/100 

0  m 

993/1000 

982/1000 

250m 

100/100 

100/100 

500m 

100/100 

100/100 

Table  1:  Example  2.  Number  of  successful  runs  (numerator)  out  of  100  or  1000  Monte  Carlo  runs 
( denominator)  for  various  values  of  y-shift  d. 


14 


Figure  3:  Example  1.  Performance  of  the  IMM/JPDAF  of  [9]  based  on  successful  (982)  runs 
(read  left  to  right,  top  to  bottom):  (a)  RMSE  in  position,  (b)  RMSE  in  velocity,  (c)  RMSE  in 
acceleration,  (d)  CV  model  probability  P\Ml{r)\Zi]  for  r  =  1,2.  (RMSE  =  root  mean-square 
error;  CV  =  constant  velocity) 


The  trajectoiy  of  the  two  targets  (y-ahift  s  -SOOm) 


The  trajectory  of  the  two  targets  (y-shift  s  50(hn) 


Figure  4:  Example  2.  The  true  trajectories  of  the  maneuvering  targets  in  terms  of  position  in 
the  xy-plane  for  two  values  of  y-shift  d.  (read  left-to-right,  top-to-bottom):  (a)  (i=-500m,  (b) 
d=500m.  Performance  (RMSE  in  position)  of  the  proposed  IMM/JPDACF  and  the  IMM/JPDAF 
of  [9]  based  on  successful  runs  for  two  values  of  y-shift  d  (read  left  to  right,  top  to  bottom):  (c) 
ci=-500m,  IMM/JPDACF,  (d)  d=500m,  IMM/JPDACF,  (e)  (i=-500m,  IMM/JPDAF,  (f)  d=500m, 
IMM/JPDAF. 


